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SUMMARY 


The major activities of the past six-month period were as follows: 

• Preparation of a paper entitled, Simulation and Analysis of a Geopotential 
Research Mission. The paper is given in Appendix A and will be submitted 
for journal publication. 

• Preparation of magnetic tapes containing the data created for the GRM 
simulation. The format of the tape is given in Appendix B. 

• Initial preparation for follow-on studies. 

On May 8, 1986, tapes containing the simulated data were delivered to Frank Lerch (GSFC) and 
Carl Wagner (NGS). On May 15, a tape was mailed to Dr. S. Boze (Applied Sciences Analytics). 

A paper entitled, Simulation of a Geopotential Research Mission (GRM), will be presented at the 
COSPAR meeting in Toulouse, France, on July 4, 1986. 
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ABSTRACT 

A computer simulation has been performed for a Geopotential Research Mission (GRM) to 
enable study of the gravitational sensitivity of the range-rate measurement between two satellites and 
to provide a set of simulated measurements to assist in the evaluation of techniques developed for 
the determination of the gravity field. The simulation, identified as SGRM 8511, was conducted 
with two satellites in near circular, frozen orbits at 160 km altitude and separated by 300 km. High 
precision numerical integration of the polar orbits was used with a gravitational field complete to 
degree and order 180 coefficients and to degree 300 in orders 0 to 10. The set of simulated data for a 
mission duration of about 32 days was generated on a Cray X-MP computer. The characteristics of 
the simulation and the nature of the results are described in this report. 

INTRODUCTION 

The determination of the earth’s gravity field by the Geopotential Research Mission (GRM) will 
be achieved using two satellites moving in essentially coincident orbit planes but separated by 100 to 
600 km. Precise measurements of the range-rate between the two polar inclination satellites would 
provide a data set from which the gravity information would be extracted. With a range-rate meas- 
urement precision of one micron per second or better and an orbit altitude of 160 km, it is expected 
that the data set would provide resolutions of at least 1-2 mgal in gravity anomalies, 5-10 cm in 
geoid height, and 100 km spatially [Taylor et al., 1984], Applications of this mission to geophysics, 
geodesy and oceanography are discussed by the National Research Council [1979] and Keating et 
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al„ [1986]. 


For the desired resolution, it is essential that the gravity field be recovered to a high degree and 
order in the spherical harmonic representation or other representations. However, use of the range- 
rate data in the procedures commonly used in satellite geodesy for gravity field determinations (e.g., 
Lerch et ah, 1985) would result in the need to solve very large systems of equations. For example, a 
degree and order 180 (or 180x180) field would require the inversion of a full, symmetric matrix of 
dimension exceeding 32,000. The computational effort associated with such an inversion would be 
significant, even for current supercomputers. Although supercomputers are expected to be available 
by 1990 with physical memory storage exceeding 1 billion 64-bit words, thereby enabling storage of 
a complete 32000x32000 matrix in main memory, the computational problems associated with such 
an inversion are still significant. 

Alternate procedures have been proposed to circumvent the need for inversion of a very large 
matrix. These procedures include those proposed by Colombo [1984], Kaula [1983], Wagner 
[1983], and others. While these procedures are promising, none have been tested on a set of simu- 
lated data derived from the solution of dynamically consistent equations of motion. 

The goal of the investigation described in this report was the creation of a set of simulated 
range-rate data that could be used for studies of techniques applicable to recovery of the gravity field 
from GRM range-rate data. The simulation model, the supporting studies of the computational accu- 
racy, and the characteristics of the resulting simulated data are discussed in the following sections. 

GRM DESCRIPTION 

The 160 km altitude proposed for the GRM results in significant forces acting on each spacecraft 
due to atmospheric resistance. Past experience has demonstrated difficulty in adequately modeling 
this force, however, its signature is very similar to that produced by various gravitational 
coefficients. As a consequence, it is essential that the force be directly removed or measured. The 
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GRM uses a "drag-free" concept, as described by Taylor et al. [1984] and Keating et ai, [1986], to 
effectively remove the effects of all nongravitational forces, including solar and earth radiation pres- 
sure as well as atmospheric forces. 

The range-rate system of the GRM operates at 91 GHz and 42 GHz to enable removal of the 
ionospheric effects on the measurement. To meet the mission objectives, it is expected that the rela- 
tive velocity variations must be measured to an accuracy of 2.5 microns per second every four 
seconds. Furthermore, it is anticipated that the range between the two spacecraft would be adjusted 
between 150 and 550 km during the six-month mission duration launched from the shuttle orbiter. 
The 2600-2800 kg spacecraft would be tracked by the Defense Mapping Agency TRANET doppler 
system. 

The GRM measurement is sensitive to gravitational effects produced by a variety of sources, 
including temporal variations. Compared to higher altitude satellites at typical geodetic altitude of 
1000 km, the influence of solid earth and ocean tide effects will be enhanced, however, the direct 
luni-solar gravitational effects will be diminished. ' 

SIMULATION DESCRIPTION 

A complete simulation of the GRM is desirable to enable the testing of techniques for determin- 
ing the gravity field from the range-rate signal in a realistic environment. However, the wide range 
of forces and kinematical effects noted in the preceding section that should be included in the simu- 
lation will complicate the studies directed toward the technique tests. 

In order to enhance the primary gravitational signal and offer the recovery techniques a set of 
simulated data that is not complicated by signals from other sources, it was concluded that the initial 
simulation should exhibit only the primary features of the actual GRM. As a consequence, the tem- 
poral variations in the gravity field and direct luni-solar forces were not included. It was also 
assumed that the spacecraft moved in an environment free of nongravitational forces; however, 
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specific modeling of the thrusting associated with the drag compensation mechanism was not 
included. Furthermore, the effects of precession, nutation, polar motion and UT1 variations were not 
included. The specific model characteristics for the simulation are shown in Table 1. 

Since the specific orbit characteristics for GRM have not been selected, the adopted orbits 
resulted from discussions within the International Association of Geodesy Special Study Group 2.83, 
chaired by Dr. R. Rummel, Delft University. Although the planned mission would use different 
separation distances over a mission duration of six months, a simulation of shorter duration with a 
single average separation of 300 km was regarded as appropriate for initial technique validation. 
The specific orbit selection was made to enable a repeating ground track over the mission duration 
with an equatorial separation between the tracks that is commensurate with the desired resolution of 
the gravity field. 

The preceding considerations resulted in the orbit characteristics given in Table 1, namely a 
simulation duration of 32 sidereal days consisting of 525 revolutions of each satellite. Furthermore, 
the adopted orbit has the characteristic that the next 525 revolutions would nearly repeat the 
groundtrack of the simulated 525 revolutions. In addition, the adopted orbits are "frozen" with peri- 
gee at 90 degrees as described by Cook [1966]. The specific initial conditions for each satellite are 
given in Table 2. 

TECHNIQUES 

The primary technique used for the creation of the simulated data was numerical integration of 
the equations of motion. Extensive studies performed on the characteristics of numerical integration 
by Lundberg [1985], as well as additional supporting studies for the GRM simulation, have demon- 
strated that truncation error resulting from the adopted numerical integration order and step size 
would not introduce significant short-period, gravity-like signals in the simulated data, at least for a 
36x36 field with an integration period of 32 days. It was not possible to perform the evaluation for 
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TABLE 1. GRM SIMULATION MODEL 


Force Model 
Gravitational 

GM = 3.986013 x 10 ,4 m 3 s -2 
a e = 6378155 m 

Static: Rapp [private communication, 1985] gravity coefficients complete to 

180x180 plus order 0 to 10 coefficients to degree 300 

Temporal: None 

Luni-Solar: None 

Nongravitational 

Drag: None 

Radiation Pressure: None 

(Note: It is assumed that the drag compensation mechanism removes all 
nongravitational forces.) 

Kinematic Models 

Earth Orientation 

Angular velocity: Constant (7.2921 1585531 * 10~ 5 rad/sec) 

Initial Greenwich mean sidereal time: 1.74731 127 radians = 100.1 135613° 

k 

Polar motion: None 
Precession/nutation: None 
Orbits 

Altitude: 160 km 
Inclination: 90 degrees 
Eccentricity: 0.00114 
Perigee: 90 degrees 
Satellite separation: 300 km 

Groundtrack repeat period: 32 sidereal days (closure <10 km) 

Duration: 32 sidereal days (2,757,250.8 s) 

the full gravity field of the simulation; however, one-day studies were performed and are described 
in the next section. Future studies will be conducted to further investigate this aspect, although it is 
not expected to produce a significant signal. 

To further aid in the control of numerical integration accuracy, the formulation of the differential 
equations of motion used an Encke-type formulation in which a reference orbit was adopted and the 
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TABLE 2. INITIAL CONDITIONS FOR SIMULATION* 


(t = 0 ) 


Lead Satellite: 
x = 0.0 

y = -150000.0 
2 =6514763.771449 
x =0.0 

y =-7817.14687496 
i =-179.5130619235 


Trailing Satellite: 

x =0.0 
y = 150000.0 
2 =6514766.990461 
x =0.0 

y =-7817.14687496 
i = 179.5130619235 


Separation distance = 300000.000017 

The initial conditions were provided by Dr. O. Colombo, EG&G Washing- 
ton Analytical Services. Units are meters for position and meters per 
second for velocity. The origin of the coordinate system coincides with 
the earth’s center of mass and the coordinate system is nonrotating. 


differential equations describe departures from the reference. Studies by Lundberg [1985] demon- 
strated that high numerical accuracy could be achieved for long-term integrations without 
rectification. In particular, the equations of motion of either satellite are 

( 1 ) 

r 3 

where F is the geocentric position vector of the satellite expressed in a nonrotating coordinate sys- 
tem, GM is the gravitational parameter of the earth, and f is the perturbing force derived from the 
gradient of the gravitational potential. The Encke procedure uses a reference orbit, F* , such that 

F = F* + e (2) 

where £ is referred to as the Encke-vector representing the satellite displacement from the reference 
orbit and described in terms of a nonrotating coordinate system. It follows that 


* 

e = -GM —z — ~TT +/-J 
r i r 3 


where s is the perturbing force on the two-body motion of the reference orbit. The advantage of this 


formulation is that by properly choosing the reference, orbit the Encke-vector would remain much 
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smaller than r (e.g., <1%), thereby resulting in more significant digits in the computation than Eq. 
(1). The numerical integration parameters used for the simulation are given in Table 3. 


TABLE 3. GRM SIMULATION INTEGRATION PARAMETERS 
Solution Method 

Encke formulation [Lundberg, 1985] 

Class 2, fixed-mesh, multistep algorithm [Lundberg, 1981] 

Step Size: 5 seconds 
Order: 10 


The accuracy of the numerical integration is dependent on the step size and order of the integra- 
tor as well as the numerical properties of the computation of the perturbing force, / . The spherical > 
harmonics formulation was adopted for the computation of the gravitational force because of its 
well-understood influence on the motion of a satellite (e.g., Kaula, 1966). The adopted formulation 
of spherical harmonics was derived by Pines [1973], a formulation with the advantage of retaining a 
close relationship with the spherical harmonic formulation while eliminating the singularities that 
exist at the poles with spherical coordinates. Because of the polar orbits used for the mission, the 
elimination of the singularities was essential. The modified Legendre functions in Pines formulation 
are computed by recursions in normalized form, hence the normalized spherical harmonic 
coefficients are used directly. Numerical experiments have been performed with possible recursions 
by computing the functions to degree and order 300 in 14 decimal digit precision and comparing 
with 28-digit results. The specific recursion algorithm adopted for the simulation is a column-wise 
formulation that has demonstrated excellent stability and numerical properties [Lundberg and 
Schutz, 1986]. The Pines formulation code has been further validated by comparison with indepen- 
dent code using the traditional formulation. 

In summary, numerical integration was used to obtain the set of simulated data. Each satellite 
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was treated independently in the process, and the simulated range and range-rate between the satel- 
lites were computed from the positions and velocities of the respective satellites. Each satellite’s 
equations of motion were formulated in an Encke-type representation to foster high-accuracy results. 

PRELIMINARY RESULTS 

Computational experiments were performed in 1983 to provide data from which an estimate of 
the computer time required for a GRM simulation could be made. This experiment consisted of 
measurements of the computer time required for the numerical evaluation of a degree and order 180 
(180x180) spherical harmonic gravity field. This early evaluation suggested that, while the formula- 
tion could be readily adapted to the vectorization environments of modem supercomputers, the time 
required for a complete simulation would be 10-30 hours; however, the time required on other main- 
frames would be 1000 hours or more. Further details on these studies with a variety of computers 
are given by Schutz [1986]. 

Preliminary simulation results for initial test and evaluation purposes were computed using a 
one-day interval of time. These results were obtained using a Cray X-MP/48 and an Amdahl 1200. 
The results from the two machines agreed to 13 decimal digits at the end of the integration period. 
The execution of the same code required 1634 seconds on the Amdahl 1200 and 1575 seconds on the 
Cray X-MP using a single processor. Furthermore, the Fortran compilers of both machines automat- 
ically vectorized the same loops in the code and both identified the recursive algorithms used in the 
computation of the modified Legendre functions as an area requiring special treatment for vectoriza- 
tion. 

The one-day preliminary results were compared with an integration performed using the simula- 
tion gravity field but truncated at degree and order 36. The initial conditions in Table 2 were used, 
and the computations with the 36x36 field were performed using a Cyber 170/750. The Cray X-MP 
result with the Table 1 model was used to compute the instantaneous range and range-rate between 
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the satellites, from which the range and range-rate computed from the 36x36 field was subtracted. 
The instantaneous range vector, p, is computed as 

P = r L - r T (4) 

and the range is p = [p • p]' A , where r T is the position vector of the trailing satellite and T L refers to 
the leading satellite. The instantaneous range-rate, p, is 



In addition to the computational speed, the precision of various machines is of importance. 
Computational precision is dependent on the number of bits allocated to storage of the floating point 
mantissa and on the central processor design which includes hardware rounding operations. The 
mantissa bit allocation for the Amdahl 1200 is 56 bits and the Cray X-MP is 48 bits. Although the 
roundoff error of the Amdahl 1200 is about 10 -15 and about 10 -14 on the Cray X-MP, effective 
hardware accuracy guards are provided on the Cray computer to assist in the control of roundoff 
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error. Experiments by the authors on various problems, including unstable dynamical systems that 
exhibit significant response to roundoff error, has shown these guards to be effective. Such accuracy 
guards were invoked in all computations using the Cray X-MP described in this report. 

SIMULATION RESULTS 

The simulation results for the full 32-day period were performed using a Cray X-MP/216 com- 
puter with a time grant from Cray Research, Incorporated (CRI). The machine was located at the 
CRI Mendota Heights, Minnesota, facility and the computations were performed in November 1985 
using a block of dedicated machine time. A subsequent verification test was performed on the Cray 
X-MP/24 installed at the University of Texas System Center for High Performance Computing in 
March 1986. The central processor time required on the Mendota Heights machine was 6 hours 28 
minutes. Furthermore, comparison with the earlier results illustrates that the vectorization of the 
Legendre recursion substantially influenced the overall computation time. Other analyses have 
shown that the gravity force evaluation consumed more than 90% of the computation time. 

The accuracy associated with the integration step size of 5 seconds was evaluated by comparison 
with an integration performed with a 2.5-second step size. However, the integration was performed 
for one day since the full 32-day computation was estimated to require about 13 hours for a step size 
of 2.5 seconds. The difference between the range and range-rate computed with the different step 
sizes are shown in Figures 4 and 5. 

The effect of step size on the range between the two satellites is illustrated in Figure 4. Over the 
one day interval, the effect is less than ten microns. More importantly, the periodic effect has an 
amplitude of about 3 microns at a frequency of one cycle per orbital revolution. Similarly, the step 
size influence on range-rate is shown in Figure 5 with a peak of less than 0.3 microns per second, but 
with a root mean square of less than 0.1 microns per second and no evidence of a linear trend in the 
differences. Consequently, it can be tentatively concluded that the 5-second step size does not 
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adversely influence the results at a level exceeding the anticipated measurement precision. 
Nevertheless, plans are being made to completely verify this conclusion through computations for 
the entire 32 days. 

The state vectors of each satellite after 32 sidereal days are given in Table 4 for the nearest 
integration node point, which differs from the exact 32 sidereal day time by 0.8 seconds. Com- 
parison with Table 2 shows that the y -components have the largest difference, about 5400 m for the 
lead satellite and 2300 m for the trailing satellite. Extrapolating the position forward in time to 
account for the 0.8-second effect results in a smaller difference for the lead satellite but a larger one 
for the trailing satellite. In either case, a closure to less than 10 km appears to have been achieved. 


TABLE 4. FINAL CONDITIONS FOR SIMULATION 


(Time = 2757250.0 seconds) 


Lead Satellite: 

x =-25.5702759 
y =-144635.9583872 
z =651517.4216461 
x =-1.4053666402 
y =-7817.7289670476 
i =-173.1909583647 


Trailing Satellite: 

x = 28.7872855 
y = 152283.2962437 
z =6514295.7534658 
x =-1.4535883592 
y =-7817.5867963785 
i = 182.1378972218 


Separation distance = 296919.342351 


Units are meters for position and meters per second for velocity. 


The time histories of the instantaneous range and range-rate for the complete 32-day simulation 
period are shown in Figures 6 and 7. The range illustrated in Figure 6 exhibits a short period (orbital 
period) with an amplitude of 450 m. Long period effects are evident also, due possibly to reso- 
nances, as well as indications that the satellites begin a linear drift at about 16 days. The range-rate 
in Figure 7 exhibits essentially a zero mean with a short period amplitude of about 0.8 m s -1 . 

Since one of the primary purposes of the GRM is to measure the gravity field, it is reasonable to 
expect that the field available at GRM launch for data analysis will be in error. Based on the planned 
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TRANET tracking system, these ground-based data would be used to determine an orbit that fit the 
data in accordance with some criteria, such as least squares. To simulate this process, the ephem- 
erides generated on the Cray X-MP with the Table 1 model were used as "observations" to an esti- 
mation process that determined epoch position and velocity of the respective satellites. In this esti- 
mation process, however, the reference gravity field differed from the Table 1 field. 

Two reference fields have been used in the investigation. In the first case, the influence of high 
degree and order coefficients was investigated by truncating the Table 1 field at degree and order 36. 
Consequently, this case used the "true" coefficients for the 36x36 portion of the field, but mismo- 
deled all coefficients of degree and order greater than 36 by effectively setting those coefficients to 
zero. In the second case, the reference field was GEM-10B [Lerch et al., 1981] which results in a 
mismodeling of all coefficients, including zero values for coefficients of degree and order exceeding 
36. In both cases, the initial position and velocity vectors (six parameters) and a constant along- 
track force parameter were estimated. The latter parameter was used to remove any inherent numeri- 
cal integration differences, however, the estimated quantity was small and approaching the level of 
computational precision indicating close compatibility between the full model integration and the 
reference model integration. 

The ephemerides resulting from the estimation using the truncated (36x36) Table 1 reference 
field were compared to the full model (Table 1) ephemerides. The position differences over 32 days 
result in position residuals with an overall RMS of 949 m for the leading satellite and 904 m for the 
trailing satellite. These differences for the leading satellite, resolved into radial, along-track and 
cross-track components, are shown in Figures 8, 9 and 10, respectively. Similar characteristics are 
exhibited by the trailing satellite. The longest period variation in the along-track component is con- 
jectured to be associated with the odd-degree zonal harmonics, while the approximately six-day 
period is believed to be due to resonant coefficients. Experiments with the second case, the GEM- 
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10B field, show position residuals with an RMS of 950 m for the leading satellite and 907 m for the 
trailing satellite, as illustrated in Figures 11, 12 and 13. Because of the general similarity between 
these two cases, it is apparent that the coefficients of degree and order greater than 36 have 
significant orbit perturbation influence on the low altitude GRM. This behavior will have some bear- 
ing on the orbit determination accuracy requirement of 300 m [Keating et al., 1986]. 

Range and range-rate residuals were computed using the Cray ephemerides resulting from the 
integration with the Table 1 model and the GEM-10B ephemerides fit to the Cray result. The range 
residuals are given in Figure 14, and the range-rate residuals are in Figure 15. Since the linear trend 
apparent in Figure 6 has either been removed or diminished, it can be concluded that the drift is pri- 
marily the result of the adopted orbit characteristics and the low degree and order coefficients. In 
general, a rich spectrum of frequencies appears to exist in Figures 14 and 15 that will provide a use- 
ful basis for further analysis, especially in the application of techniques designed to recover the grav- 
ity field from such data. 

RECOMMENDATIONS 

This report has described the results of the first simulation of a GRM with emphasis on the gravi- 
tational aspects. The simulation model has been chosen to enhance the gravitational signal for the 
testing and validation of gravity recovery techniques. The results of the simulation are available for 
further analysis on a high density (6250 cpi) magnetic tape. 

Based on the results described in this report, several recommendations for future study can be 
made as follows: 

• The orbit determination accuracy will be significantly influenced by a priori errors in the gravity 
field, including coefficients exceeding degree and order 36; consequently techniques should be 
investigated that reduce the influence of the resulting orbit errors on the GRM gravity field 
recovery. . 
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• A more realistic simulation of the TRANET tracking distribution should be investigated directed 
toward the influence on the range-rate residuals. 

• Examination of alternate tracking systems, such as satellite-borne Global Positioning System 
receivers should be performed. 

• Implementation of more realistic force models, especially tide models, and kinematic models, 
especially those describing the orientation of the earth in space. 

• Simulation of outer satellite drag and creation of a set of range and range-rate simulated data that 
include thrusts associated with a drag compensation mechanism. 

• Investigation of the influence of time-varying features in the ocean surface on the range and 
range-rate. 

• Further investigation of the numerical precision associated with the simulation. 

• Examination of the influence of satellite separation and altitude. 

Efforts to investigate these GRM aspects are in progress. 
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Figure 4. Range residuals, 5-second minus 2.5-second step sizes, full model for one day 
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Figure 5. 


Range-rate residuals, 5-second minus 2.5-second step sizes, full model for one day 
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Figure 6. Range history for Table 1 model 
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Figure 7. Range-rate history for Table 1 model 
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APPENDIX B 


DISTRIBUTION FORMAT 



Center for Space Research 
The University of Texas at Austin 
May 2, 1986 

GRM SIMULATION: SGRM8511 
Distribution Format 


Characters per record: 80 

Records per block: 300 

Characters per block: 24000 

Character set: either ASCII or EBCDIC (as noted on label) 

Density: 6250 cpi 

There are three files on the tape: 

File #1: Contains the time, range, and range-rate from the GRM simulation using the OSU322 

gravity (180x180 and orders 0 to 10 for degrees 181 to 300). Also contains the 
ephemerides for a reference orbit fit to the GRM simulation using the modified GEM- 
10B geopotential. 

File #2: OSU322 gravity file. 

File #3: GEM-10B gravity file, with modified GM and a e 

The formats for reading files 2 and 3 are at the beginning for each file. The format and data content 
for the first file is as follows: 

Data at each time point: 

f,- time tag in seconds from the start of the ephemeris at t = 0. 
pi observed range between the two satellites (noiseless) 
pi observed range-rate between the two satellites (noiseless) 
r|i * measurement noise parameter 

7 ] ** position vector of the lead satellite for the 36x36 reference orbit 
v i ** velocity vector of the lead satellite for the 36x36 reference orbit 
7 2 ** position vector of the trailing satellite for the 36x36 reference orbit 
v 2 ** velocity vector of the trailing satellite for the 36x36 reference orbit 

* The measurement noise parameter ranges from -5.0 to +5.0 and has a Gaussian distribution. 
Scaling this quantity by the standard error of the measurement precision enables noise to be 
added to the observations. For example, to create a range-rate observation with noise, p,-, 
p,- = p,- + or),-, where a might be 0.000001 m s _1 . . 
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** These data are expressed in an earth-fixed, rotating coordinate system. The position and velocity 
vectors were transformed from the nonrotating, inertial coordinate system of integration to the 
earth-fixed system using the same transformation matrix, i.e., 

n = [R]r fi 
V; = 

where 7[ t and v}. are the position and velocity vectors in inertial coordinates and [7?] is the 
transformation matrix between the nonrotating, inertial system and the earth-fixed, rotating 
system. 

Three records (card images) are used to display the data at each time point; the card images and 
formats are: 
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The ground tracks of the satellites nearly repeat after 32 sidereal days (2757250.8 sec); the data 
stored on the magnetic tape span 32 mean solar days (2764800 sec). 
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